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Abstract 

In  this  paper,  a pore-scale  network  modeling  method,  based  on  the  flow  con- 
tinuity residual  in  conjunction  with  a Newton-Raphson  non-linear  iterative 
solving  technique,  is  proposed  and  used  to  obtain  the  pressure  and  flow  fields 
in  a network  of  interconnected  distensible  ducts  representing,  for  instance, 
blood  vasculature  or  deformable  porous  media.  A previously  derived  an- 
alytical expression  correlating  boundary  pressures  to  volumetric  flow  rate 
in  compliant  tubes  for  a pressure-area  constitutive  elastic  relation  has  been 
used  to  represent  the  underlying  flow  model.  Comparison  to  a preceding 
equivalent  method,  the  one-dimensional  Navier-Stokes  finite  element,  was 
made  and  the  results  were  analyzed.  The  advantages  of  the  new  method 
have  been  highlighted  and  practical  computational  issues,  related  mainly  to 
the  rate  and  speed  of  convergence,  have  been  discussed. 

Keywords:  fluid  mechanics;  one-dimensional  flow;  Navier-Stokes;  distensible 
network;  compliant  porous  media;  blood  circulation;  non-linear  system. 


1 Introduction 

There  are  many  scientific,  industrial  and  biomedical  applications  related  to  the 
flow  of  fluids  in  distensible  networks  of  interconnected  tubes  and  compliant  porous 
materials.  A few  examples  are  magma  migration,  microfluidic  sensors,  fluid  filtering 
devices,  deformable  porous  geological  structures  such  as  those  found  in  petroleum 
reservoirs  and  aquifers,  as  well  as  almost  all  the  biological  flow  phenomena  like 
blood  circulation  in  the  arterial  and  venous  vascular  trees  or  biological  porous 
tissue  and  air  movement  in  the  lung  windpipes. 

There  have  been  many  studies  in  the  past  related  to  this  subject  [1-10];  however 
most  of  these  studies  are  based  on  complex  numerical  techniques  built  on  tortuous 
mathematical  infrastructures  which  are  not  only  difficult  to  implement  with  expen- 
sive computational  running  costs,  but  are  also  difficult  to  verify  and  validate.  The 
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widespread  approach  in  modeling  the  flow  in  deformable  structures  is  to  use  the 
one- dimensional  Navier-Stokes  finite  element  formulation  for  modeling  the  flow  in 
networks  of  compliant  large  tubes  [5,  11]  and  the  extended  Darcy  formulation  for 
the  flow  in  deformable  porous  media  which  is  based  on  the  poromechanics  theory 
or  some  similar  numerical  meshing  techniques  [12-15].  Rigid  network  flow  models, 
like  Poiseuille,  may  also  be  used  as  an  approximation  although  in  most  cases  this 
is  not  really  a good  one  [16]. 

There  have  also  been  many  studies  in  the  past  related  to  the  flow  of  fluids  in 
ensembles  of  interconnected  ducts  using  pore-scale  network  modeling  especially  in 
the  earth  science  and  petroleum  engineering  disciplines  [17-23].  However,  there  is 
hardly  any  work  on  the  use  of  pore-scale  network  modeling  to  simulate  the  flow 
of  fluids  in  deformable  structures  with  distensible  characteristics  such  as  elastic  or 
viscoelastic  mechanical  properties. 

There  are  several  major  advantages  in  using  pore-scale  network  modeling  over 
the  more  traditional  analytical  and  numerical  approaches.  These  advantages  in- 
clude a comparative  ease  of  implementation,  relatively  low  computational  cost, 
reliability,  robustness,  relatively  smooth  convergence,  ease  of  verification  and  val- 
idation, and  obtaining  results  which  are  usually  very  close  to  the  underlying  ana- 
lytical model  that  describes  the  flow  in  the  individual  ducts.  Added  to  all  these  a 
fair  representation  and  realistic  description  of  the  flow  medium  and  the  essential 
physics  at  macroscopic  and  mesoscopic  levels  [24,  25].  Pore-scale  modeling,  in  fact, 
is  a balanced  compromise  between  the  technical  complexities  and  the  physical  re- 
ality. More  details  about  pore-scale  network  modeling  approach  can  be  found,  for 
instance,  in  [26]. 

In  this  paper  we  use  a residual-based  non-linear  solution  method  in  conjunction 
with  an  analytical  expression  derived  recently  [27]  for  the  one- dimensional  Navier- 
Stokes  flow  in  elastic  tubes  to  obtain  the  pressure  and  flow  fields  in  networks  of 
interconnected  distensible  ducts.  The  residual-based  scheme  is  a standard  method 
for  solving  systems  of  non-linear  equations  and  hence  is  commonly  used  in  fluid  me- 
chanics for  solving  systems  of  partial  differential  equations  obtained,  for  example, 
in  a finite  element  formulation  [11],  The  proposed  method  is  based  on  minimizing 
the  residual  obtained  from  the  conservation  of  volumetric  flow  rate  on  the  indi- 
vidual network  nodes  with  a Newton-Raphson  non-linear  iterative  solution  scheme 
in  conjunction  with  the  aforementioned  analytical  expression.  Other  analytical, 
empirical  and  even  numerical  relations  [28]  describing  the  flow  in  deformable  ducts 
can  also  be  used  to  characterize  the  underlying  flow  model. 
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2 Method 


The  flow  of  an  incompressible  Newtonian  fluid  in  a tube  with  length  L and  cross 
sectional  area  A assuming  a laminar  axi-symmetric  slip-free  flow  with  a fixed  profile 
and  negligible  gravitational  body  forces  can  be  described  by  the  following  one- 
dimensional Navier-Stokes  system  of  mass  and  momentum  conservation  relations 


where  Q is  the  volumetric  flow  rate,  t is  the  time,  x is  the  axial  coordinate  along 
the  tube,  a is  the  momentum  flux  correction  factor,  p is  the  fluid  mass  density,  p 
is  the  local  pressure,  and  n is  the  viscosity  friction  coefficient  which  is  normally 
given  by  k = with  v being  the  fluid  kinematic  viscosity  defined  as  the  ratio 
of  the  dynamic  viscosity  p to  the  mass  density  [11,  16,  29-31].  These  relations 
are  usually  supported  by  a constitutive  relation  that  correlates  the  pressure  to  the 
cross  sectional  area  in  a distensible  tube,  to  close  the  system  in  the  three  variables 
A,  Q and  p. 

The  usual  method  for  solving  this  system  of  equations  for  a single  compliant 
tube  in  transient  and  steady  state  flow  is  to  use  the  finite  element  method  based 
on  the  weak  formulation  by  multiplying  the  mass  and  momentum  conservation 
equations  by  weight  functions  and  integrating  over  the  solution  domain  to  obtain 
the  weak  form  of  the  system.  This  weak  form,  with  suitable  boundary  conditions, 
can  then  be  used  as  a basis  for  finite  element  implementation  in  conjunction  with 
an  iterative  scheme  such  as  Newton-Raphson  method.  The  finite  element  system 
can  also  be  extended  to  a network  of  interconnected  deformable  tubes  by  imposing 
suitable  boundary  conditions,  based  on  pressure  or  flux  constraints  for  instance, 
on  all  the  boundary  nodes,  and  coupling  conditions  on  all  the  internal  nodes.  The 
latter  conditions  are  normally  derived  from  Riemann’s  method  of  characteristics, 
and  the  conservation  principles  of  mass  and  mechanical  energy  in  the  form  of 
Bernoulli  equation  for  inviscid  flow  with  negligible  gravitational  body  forces  [32], 
More  details  on  the  finite  element  formulation,  validation  and  implementation  are 
given  in  [11]. 

The  pore-scale  network  modeling  method,  which  is  proposed  as  a substitute  for 
the  finite  element  method  in  steady  state  flow,  is  established  on  three  principles: 
the  continuity  of  mass  represented  by  the  conservation  of  volumetric  flow  rate  for 
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incompressible  flow,  the  continuity  of  pressure  where  each  branching  nodal  point 
has  a uniquely  defined  pressure  value  [32],  and  the  characteristic  relation  for  the 
flow  of  the  specific  fluid  model  in  the  particular  structural  geometry  such  as  the 
flow  of  power  law  fluids  in  rigid  tubes  or  the  flow  of  Newtonian  fluids  in  elastic 
ducts.  The  latter  principle  is  essentially  a fluid-structure  interaction  attribute  of 
the  adopted  flow  model  especially  in  the  context  of  compliant  ducts. 

In  more  technical  terms,  the  pore-scale  network  modeling  method  employs  an 
iterative  scheme  for  solving  the  following  matrix  equation  which  is  based  on  the 
flow  continuity  residual 


JAp  = — r (3) 

where  J is  the  Jacobian  matrix,  p is  the  vector  of  variables  which  represent  the 
pressure  values  at  the  boundary  and  branching  nodes,  and  r is  the  vector  of  resid- 
uals which  is  based  on  the  continuity  of  the  volumetric  flow  rate.  For  a network  of 
interconnected  tubes  defined  by  n boundary  and  branching  nodes  the  above  matrix 
equation  is  defined  by 
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where  the  subscripts  stand  for  the  nodal  indices,  p and  r are  the  nodal  pressure 
and  residual  respectively,  and  / is  the  flow  continuity  residual  function  which,  for 
a general  node  j,  is  given  by 


f,  = J2Qi  = ° (5) 

i= 1 

In  the  last  equation,  m is  the  number  of  flow  ducts  connected  to  node  j,  and  Qi 
is  the  volumetric  flow  rate  in  duct  i signed  (+/— ) according  to  its  direction  with 
respect  to  the  node,  i.e.  toward  or  away.  For  the  boundary  nodes,  the  continuity 
residual  equations  are  replaced  by  the  boundary  conditions  which  are  usually  based 
on  the  pressure  or  flow  rate  constraints.  In  the  computational  implementation,  the 
Jacobian  is  normally  evaluated  numerically  by  finite  differencing  [11]. 

The  procedure  to  obtain  a solution  by  the  residual-based  pore-scale  modeling 
method  starts  by  initializing  the  pressure  vector  p with  initial  values.  Like  any 
other  numerical  technique,  the  rate  and  speed  of  convergence  is  highly  dependent 
on  the  initial  values  of  the  variable  vector.  The  system  given  by  Equation  4 is  then 
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constructed  where  the  Jacobian  matrix  and  the  residual  vector  are  calculated  in 
each  iteration.  The  system  4 is  then  solved  for  Ap,  i.e. 

Ap  = — J_1r  (6) 

and  the  vector  p in  iteration  l is  updated  to  obtain  a new  pressure  vector  for  the 
next  iteration  (7  + 1),  that  is 


Pm  = pi  + Ap  (7) 

This  is  followed  by  computing  the  norm  of  the  residual  vector  from  the  following 
equation 


m = + • • • + r2n 


(8) 


n 

where  r is  the  flow  continuity  residual.  This  cycle  is  repeated  until  the  norm  is  less 
than  a predefined  error  tolerance  or  a certain  number  of  iteration  cycles  is  reached 
without  convergence.  In  the  last  case,  the  operation  will  be  deemed  a failure  and 
hence  it  will  be  aborted  to  be  resumed  possibly  with  improved  initial  values  or 
even  modified  model  parameters  if  the  physical  problem  is  flexible  and  allows  for 
a certain  degree  of  freedom. 

The  characteristic  flow  relation  that  has  to  be  used  for  computing  Q in  the  resid- 
ual equation  is  dependent  on  the  flow  model.  As  for  the  flow  of  Newtonian  fluids  in 
distensible  tubes  based  on  the  previously-described  system  of  flow  equations,  the 
following  analytical  relation  representing  the  one-dimensional  Navier-Stokes  flow 
in  clastic  tubes  can  be  used 


-x-L  + J KW  - 4a  In  (Am/Am)  +2) 

o = Y 2 A (O') 

2a  In  (Ain/Aou) 

Other  analytical  or  empirical  or  numerical  relations  characterizing  the  flow  rate 
can  also  be  used  in  this  context  [28]. 

The  flow  relation  of  Equation  9 was  previously  derived  and  validated  by  a one- 
dimensional finite  element  method  in  [27].  Equation  9 is  based  on  a pressure-area 
constitutive  elastic  relation  in  which  the  pressure  is  proportional  to  the  radius 
change  with  a proportionality  stiffness  factor  that  is  scaled  by  the  reference  area, 

i.e. 
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In  the  last  two  equations,  Aa  is  the  reference  area  corresponding  to  the  reference 
pressure  which  in  this  equation  is  set  to  zero  for  convenience  without  affecting  the 
generality  of  the  results,  Ain  and  Aou  are  the  tube  cross  sectional  area  at  the  inlet 
and  outlet  respectively,  A is  the  tube  cross  sectional  area  at  the  actual  pressure, 
p,  as  opposed  to  the  reference  pressure,  and  f3  is  the  tube  wall  stiffness  coefficient 
which  is  usually  defined  by 

_ y/nhoE 

P 1 9 

l - <r 

where  hQ  is  the  tube  wall  thickness  at  reference  pressure,  while  E and  ^ are  respec- 
tively the  Young’s  elastic  modulus  and  Poisson’s  ratio  of  the  tube  wall. 

With  regard  to  the  validation  of  the  numeric  solutions  obtained  from  the  fi- 
nite element  and  pore-scale  methods,  the  time  independent  solutions  of  the  one- 
dimensional finite  element  model  can  be  tested  for  validation  by  satisfying  the 
boundary  and  coupling  conditions  as  well  as  the  analytic  solution  given  by  Equa- 
tion 9 on  each  individual  duct,  while  the  solutions  of  the  residual-based  pore-scale 
modeling  method  are  validated  by  testing  the  boundary  conditions  and  the  conti- 
nuity of  volumetric  flow  rate  at  each  internal  node,  as  well  as  the  analytic  solution 
given  by  Equation  9 which  is  inevitably  satisfied  if  the  continuity  equation  is  sat- 
isfied according  to  the  pore-scale  solution  scheme.  The  necessity  to  satisfy  the 
analytic  solution  on  each  individual  tube  in  the  finite  element  method  is  based 
on  the  fact  that  the  flow  in  the  individual  tubes  according  to  the  underlying  one- 
dimensional model  is  dependent  on  the  imposed  boundary  conditions  but  not  on 
the  mechanism  by  which  these  conditions  are  imposed.  In  the  case  of  finite  el- 
ement with  tube  discretization  and/or  employing  non-linear  interpolation  orders, 
the  solution  at  the  internal  points  of  the  ducts  can  also  be  tested  by  satisfying  the 
following  analytical  relation  [11] 

-aQ2  In (A/A-m)  + /3  (a5/2  - Afn/2)  /(5 pAa) 

X = -[2irais/(a-l)\Q  (12) 

The  derivation  of  this  equation  is  similar  to  the  derivation  of  Equation  9 but  with 
using  the  inlet  boundary  condition  only.  In  fact  even  Equation  9 can  be  used  for 
testing  the  solution  at  the  internal  points  if  we  assume  these  points  as  periphery 
nodes  [11]. 
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3 Implementation  and  Results 


The  residual-based  pore-scale  modeling  method,  as  described  in  the  last  section, 
was  implemented  in  a computer  code  with  an  iterative  Newton-Raphson  method 
that  includes  four  numeric  solvers  (SPARSE,  SUPERLU,  UMFPACK,  and  LA- 
PACK).  The  code  was  then  tested  on  computer-generated  networks  representing 
distensible  fluid  transportation  structures  like  ensembles  of  interconnected  tubes 
or  porous  media.  A sample  of  these  networks  are  given  in  Figure  1.  Because  the 
residual-based  pore-scale  method  can  be  used  in  general  to  obtain  flow  solutions  for 
any  characteristic  flow  that  involves  linear  or  non-linear  fluid  models,  such  as  New- 
tonian or  non-Newtonian  fluids,  passing  through  rigid  or  distensible  networks,  the 
code  was  tested  first  on  Poiseuillc  and  power  law  fluids  in  rigid  networks  [16,  33]. 
The  results  for  these  validation  tests  were  exceptionally  accurate  with  very  low 
error  margin  over  the  whole  network  and  with  smooth  convergence. 

We  also  used  the  one-dimensional  finite  element  model  that  we  briefly  described 
in  the  last  section  for  the  purpose  of  comparison.  This  model  was  previously  im- 
plemented in  a computer  code  with  a residual-based  Newton-Raphson  iterative  so- 
lution scheme,  similar  to  the  one  used  in  the  pore-scale  modeling.  Full  description 
of  the  finite  element  method,  code  and  techniques  can  be  found  in  [11].  A number 
of  pore-scale  and  one-dimensional  finite  element  time  independent  flow  simulations 
were  carried  out  on  our  computer-generated  networks  using  a range  of  physical  pa- 
rameters defining  the  fluid  and  structure  as  well  as  different  numeric  solvers  with 
different  solving  schemes.  Various  types  of  pressure  and  flow  boundary  conditions 
were  imposed  in  these  simulations,  although  in  most  cases  Dirichlet-type  pressure 
boundary  conditions  were  applied.  The  finite  element  simulations  were  performed 
using  a linear  Lagrange  interpolation  scheme  with  no  tube  discretization  to  closely 
match  the  pore-scale  modeling  approach.  All  the  results  of  the  reported  and  un- 
reported runs  have  passed  the  rigorous  validation  tests  that  we  stated  earlier  in 
section  2. 

Regarding  the  nature  of  the  networks,  several  types  of  networks  have  been 
generated  and  used  in  the  above-mentioned  flow  simulations;  these  include  frac- 
tal, cubic  and  orthorhombic  networks.  The  fractal  networks  are  based  on  fractal 
branching  patterns  where  each  generation  of  the  branching  tubes  in  the  network 
have  a specific  number  of  branches  related  to  the  number  of  branches  in  the  par- 
ent generation,  such  as  2:1,  as  well  as  specific  branching  angle,  radius  branching 
ratio  and  length  to  radius  ratio.  The  radius  branching  ratio  is  normally  based  on 
a Murray-type  relation  [32,  34-36].  The  fractal  networks  are  also  characterized  by 
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the  number  of  generations.  The  fractal  networks  used  in  this  study  have  a single 
inlet  boundary  node  and  multiple  outlet  boundary  nodes  [16]. 

The  cubic  and  orthorhombic  networks  are  based  on  a cubic  or  orthorhombic 
three-dimensional  lattice  structure  where  the  radii  of  the  tubes  in  the  network 
can  be  constant  or  subject  to  statistical  random  distributions  such  as  uniform 
distribution.  These  networks  have  a number  of  inlet  boundary  nodes  on  one  side 
and  a similar  number  of  outlet  boundary  nodes  on  the  opposite  side  while  the 
nodes  on  the  other  four  sides  (i.e.  the  lateral)  are  considered  internal  nodes.  The 
boundary  conditions  are  then  imposed  on  these  inlet  and  outlet  boundary  nodes 
individually  according  to  the  need.  A sample  of  the  fractal,  cubic  and  orthorhombic 
networks  used  in  this  investigation  are  shown  in  Figure  1. 

ft  is  noteworthy  that  all  the  networks  used  in  the  pore-scale  and  finite  ele- 
ment flow  simulations  consist  of  interconnected  straight  cylindrical  tubes  where 
each  tube  is  characterized  by  a constant  radius  over  its  entire  length  and  spa- 
tially identified  by  two  end  nodes.  These  networks  are  totally  connected,  that  is 
any  node  in  the  network  can  be  reached  from  any  other  node  by  moving  entirely 
inside  the  network  space.  As  for  the  physical  size,  we  used  different  sizes  to  repre- 
sent different  flow  structures  such  as  arterial  and  venous  blood  vascular  trees  and 
spongy  porous  geological  media.  The  physical  parameters  used  in  our  simulations, 
especially  those  related  to  the  fluids  and  flow  structures,  were  generally  selected 
to  represent  realistic  physical  systems  although  physical  parameters  representing 
hypothetical  conditions  have  also  been  used  for  the  purpose  of  test  and  valida- 
tion. However,  since  the  current  study  is  purely  theoretical  with  no  involvement 
of  experimental  or  observational  data,  the  validity  of  the  reported  models  are  not 
affected  by  the  actual  values  of  the  physical  parameters  although  this  has  some 
consequences  on  the  speed  and  rate  of  convergence  in  different  physical  regimes. 

In  Figures  2 and  3 we  present  a sample  of  the  above  mentioned  comparative 
simulations.  In  Figure  2 we  plot  the  ratio  of  the  flow  rate  obtained  from  the  finite 
element  model  to  that  obtained  from  the  pore-scale  model  for  a fractal  network 
with  an  area-preserving  branching  index  of  2 [32]  where  the  number  of  tubes  in  each 
generation  is  twice  the  number  in  the  parent  generation.  In  these  flow  simulations 
we  applied  an  inlet  boundary  pressure  of  2000  Pa  on  the  single  inlet  boundary  node 
of  the  main  branch  and  an  outlet  boundary  pressure  of  0.0  Pa  on  all  the  outlet 
boundary  nodes.  The  results  of  the  pore-scale  and  finite  element  models  are  very 
close  although  the  two  models  differ  due  to  the  use  of  different  branching  coupling 
conditions,  i.e.  continuity  of  pressure  for  the  pore-scale  model  and  Bernoulli  for  the 
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finite  element.  The  effect  of  the  coupling  conditions  on  the  different  generations  of 
the  fractal  network  can  be  seen  in  Figure  2 where  a generation-based  configuration 
is  obvious.  This  feature  is  a clear  indication  of  the  effect  of  the  coupling  conditions 
on  the  deviation  between  the  two  models. 

In  Figure  3 we  compare  the  pore-scale  and  finite  element  models  for  an  inho- 
mogeneous orthorhombic  network  consisting  of  about  11000  interconnected  tubes, 
similar  to  the  one  depicted  in  Figure  1 (f),  where  we  plotted  the  ratio  of  the  flow 
rate  obtained  from  the  two  models  as  for  the  fractal  network.  The  network  is  gener- 
ated with  a random  uniform  distribution  for  the  tubes  radii  with  a variable  length 
in  different  orientations  and  with  a variable  length  to  radius  ratio  that  ranges  be- 
tween about  5-15.  The  dimensions  of  the  flow  structure  are  2 x 1.5  x 1 m with 
a constant  inlet  boundary  pressure  of  3000  Pa  applied  to  all  the  nodes  on  the  in- 
let side  and  zero  outlet  boundary  pressure  applied  to  all  the  nodes  on  the  outlet 
side.  The  fluid  and  structure  physical  parameters  for  the  flow  model  are  selected 
to  roughly  resemble  the  flow  of  crude  oil  in  some  elastic  structure,  possibly  in  a 
refinement  processing  plant.  As  seen,  the  pore-scale  and  finite  element  results  differ 
significantly  on  most  part  of  the  network.  The  reason  in  our  judgment  is  the  effect 
of  the  branching  coupling  conditions  (i.e.  continuity  of  pressure  and  Bernoulli) 
which  have  a stronger  impact  in  such  a network  than  a fractal  network  due  to  the 
inhomogeneity  with  random  radius  distribution  on  one  hand  and  the  high  nodal 
connectivity  of  the  orthorhombic  lattice  on  the  other  hand.  The  fluid  property 
which  significantly  differs  from  that  of  the  fractal  network  simulation  may  also 
have  a role  in  exacerbating  the  discrepancy. 

The  results  shown  in  these  figures  represent  a sample  of  our  simulations  which 
reflect  the  general  trend  in  other  simulations.  However  the  agreement  between  the 
pore-scale  and  finite  element  models  is  highly  dependent  on  the  flow  regime  and 
the  nature  of  the  physical  problem  which  combines  the  fluid,  structure  and  their 
interaction.  As  indicated  early,  the  discrepancy  between  the  two  models  reflects 
the  effect  of  the  coupling  conditions,  i.e.  the  pressure  continuity  for  the  pore- 
scale  model  and  the  Bernoulli  condition  for  the  finite  element.  The  gravity  of  this 
effect  is  strongly  dependent  on  the  type  of  fluid,  flow  regime,  inhomogeneity  and 
connectivity. 

It  should  be  remarked  that  in  these  figures  (i.e.  2 and  3)  we  used  the  volumetric 
flow  rate,  rather  than  the  nodal  pressure,  to  make  the  comparison.  The  reason  is 
that  comparing  the  pressure  is  not  possible  because  nodal  pressure  is  not  defined 
in  the  finite  element  model  due  to  the  use  of  the  Bernoulli  equation  [32]  where  each 
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node  has  a number  of  pressure  values  matching  the  number  of  the  connected  tubes. 


4 Pore-Scale  vs.  Finite  Element 

It  is  difficult  to  make  an  entirely  fair  comparison  between  the  pore-scale  and  finite 
element  methods  due  mainly  to  the  use  of  different  coupling  conditions  at  the 
branching  junctions  as  well  as  different  theoretical  assumptions.  Therefore,  the 
pressure  and  flow  rate  fields  obtained  from  these  two  methods  on  a given  network 
are  generally  different.  The  difference,  however,  is  highly  dependent  on  the  nature 
of  the  specified  physical  and  computational  conditions. 

One  of  the  advantages  of  the  pore-scale  modeling  method  over  the  finite  element 
method,  in  addition  to  the  general  advantages  of  the  pore-scale  modeling  approach 
which  were  outlined  earlier,  is  that  when  pore-scale  method  converges  it  usually 
converges  to  the  underlying  analytic  solution  with  negligible  marginal  errors  over 
the  whole  network,  while  the  finite  element  method  normally  converges  with  signifi- 
cant errors  over  some  of  the  network  ducts  especially  those  with  eccentric  geometric 
characteristics  such  as  very  low  length  to  radius  ratio  [11],  It  may  also  be  argued 
that  the  coupling  condition  used  in  the  pore-scale  modeling  method,  which  is  based 
on  the  continuity  of  pressure,  is  better  than  the  corresponding  coupling  condition 
used  in  the  finite  element  method  which  is  based  on  the  Bernoulli  inviscid  flow 
with  discontinuous  pressure  at  the  nodal  points.  Some  of  the  criticism  to  the  use 
of  Bernoulli  as  a coupling  condition  is  outlined  in  [32] . Another  advantage  of  the 
pore-scale  modeling  is  that  it  is  generally  more  stable  than  the  finite  element  with 
a better  convergence  behavior  due  partly  to  the  simpler  pore-scale  computational 
infrastructure. 

The  main  advantage  of  the  finite  element  method  over  the  pore-scale  modeling 
method  is  that  it  accommodates  time  dependent  flow  naturally,  as  well  as  time 
independent  flow,  while  pore-scale  modeling  in  its  current  formulation  is  capable 
only  of  dealing  with  time  independent  flow.  Moreover,  the  finite  element  method 
may  be  better  suited  for  describing  other  one-dimensional  transportation  phenom- 
ena such  as  wave  propagation  and  reflection  in  deformable  networks.  However, 
time  dependent  flow  can  be  simulated  within  the  pore-scale  modeling  framework 
as  a series  of  time  independent  frames  although  this  is  not  really  a time  dependent 
flow  but  rather  a pseudo  time  dependent.  Another  advantage  of  the  finite  element 
method  is  that  it  is  capable,  through  the  use  of  segment  discretization  and  higher 
orders  of  interpolation,  of  computing  the  pressure  and  flow  rate  fields  at  the  inter- 
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(c)  homogeneous  cubic  (d)  inhomogeneous  cubic 


(e)  homogeneous  orthorhombic  (f)  inhomogeneous  orthorhombic 


Figure  1:  A sample  of  computer-generated  fractal,  cubic  and  orthorhombic  net- 
works used  in  the  current  investigation. 
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Figure  2:  The  ratio  of  flow  rate  of  finite  element  to  pore-scale  models  versus  tube 
index  for  a fractal  network.  The  network  consists  of  10  generations  with  an  area- 
preserving branching  index  of  2 [32]  and  an  inlet  main  branch  with  R — 5 mm 
and  L = 50  mm.  The  parameters  for  these  simulations  are:  f3  = 236.3  Pa.m, 
p = 1060.0  kg.m-3,  /i  = 0.0035  Pa.s,  and  a = 1.333.  The  inlet  and  outlet  pressures 
are:  pt  = 2000  Pa  and  pQ  = 0.0  Pa.  The  fluid  and  structure  parameters  are  chosen 
to  roughly  resemble  blood  circulation  in  large  vessels. 


nal  points  along  the  tubes  length  and  not  only  on  the  tubes  periphery  points  at 
the  nodal  junctions.  However,  due  to  the  incompressibility  of  the  flow,  computing 
the  flow  rate  at  the  internal  points  is  redundant  as  it  is  identical  to  the  flow  rate  at 
the  end  points.  With  regard  to  computing  the  pressure  field  on  the  internal  points, 
it  can  also  be  obtained  by  pore-scale  modeling  method  through  the  application  of 
Equation  12  to  the  solution  obtained  on  the  individual  ducts.  Moreover,  it  can  be 
obtained  by  creating  internal  nodal  junctions  along  the  tubes  through  the  use  of 
tube  discretization,  similar  to  the  discretization  in  the  finite  element  method. 

With  regard  to  the  size  of  the  problem,  which  directly  influences  the  ensuing 
memory  cost  as  well  as  the  CPU  time,  the  number  of  degrees  of  freedom  for  the 
pore-scale  model  is  half  the  number  of  degrees  of  freedom  for  the  one-dimensional 
finite  element  model  due  to  the  fact  that  the  former  has  one  variable  only  (p) 
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Figure  3:  The  ratio  of  flow  rate  of  finite  element  to  pore-scale  models  versus 
tube  index  for  an  inhomogeneous  orthorhombic  network.  The  network  consists  of 
about  11000  tubes  with  different  lengths  and  length  to  radius  ratios  as  explained 
in  the  main  text.  The  parameters  for  these  simulations  are:  f3  = 236.3  Pa.m, 
p = 860.0  kg.m-3,  (i  — 0.075  Pa.s,  and  a = 1.2.  The  inlet  and  outlet  pressures  are: 
Pi  = 3000  Pa  and  pQ  = 0.0  Pa.  The  fluid  and  structure  parameters  are  chosen  to 
approximately  match  the  transport  of  crude  oil  through  an  clastic  structure. 


while  the  latter  has  two  variables  (p  and  Q).  This  estimation  of  the  finite  element 
computational  cost  is  based  on  using  a linear  interpolation  scheme  with  no  tube 
discretization;  and  hence  this  cost  will  substantially  increase  with  the  use  of  dis- 
cretization and/or  higher  orders  of  interpolation.  The  computational  cost  for  both 
models  also  depends  on  the  type  of  the  solver  used  such  as  being  sparse  or  dense, 
and  direct  or  iterative,  as  well  as  some  problem-specific  implementation  overheads. 

CPU  processing  time  depends  on  several  factors  such  as  the  size  and  type  of  the 
network,  the  initial  values  for  the  flow  solutions,  the  parameters  of  the  fluid  and 
tubes,  the  employed  numerical  solver,  and  the  assumed  pressure-area  constitutive 
relation.  Typical  processing  time  for  a single  run  of  the  pore-scale  network  model 
on  a typical  laptop  or  desktop  computer  ranges  between  a few  seconds  to  few 
minutes  using  a single  processor  with  an  average  speed  of  2-3  gigahertz.  The 
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CPU  processing  time  for  the  time  independent  finite  element  model  is  comparable 
to  the  processing  time  of  the  pore-scale  network  model.  In  both  cases,  the  final 
convergence  in  a typical  problem  is  normally  reached  within  3-7  Newton-Raphson 
iterations  depending  mainly  on  the  initial  values. 

5 Convergence  Issues 

Like  the  one-dimensional  Navier-Stokes  finite  element  model,  the  residual-based 
pore-scale  method  may  suffer  from  convergence  difficulties  due  to  the  highly  non- 
linear nature  of  the  flow  model.  The  nonlinearity  increases,  and  hence  the  conver- 
gence difficulties  aggravate,  with  increasing  the  pressure  gradient  across  the  flow 
domain.  The  nonlinearity  also  increases  with  eccentric  values  representing  the  fluid 
and  structure  parameters  such  as  the  fluid  viscosity  or  wall  distensibility.  Several 
numerical  tricks  and  stabilization  techniques  can  be  used  to  improve  the  rate  and 
speed  of  convergence.  These  include  non-dimensionalization  of  the  flow  equations, 
using  a variety  of  unit  systems  such  as  m.kg.s  or  mm.g.s  or  m.g.s  for  the  input 
data  and  parameters,  and  scaling  the  network  flow  model  up  or  down  to  obtain  a 
similarity  solution  that  can  be  scaled  back  to  obtain  the  final  solution.  The  error 
tolerance  for  the  convergence  criterion  which  is  based  on  the  residual  norm  may 
also  be  increased  to  enhance  the  rate  and  speed  of  convergence.  Despite  the  fact 
that  the  use  of  relatively  large  error  tolerance  can  cause  a convergence  to  a wrong 
solution  or  to  a solution  with  large  errors,  the  solution  can  always  be  tested  by  the 
above-mentioned  validation  metrics  and  hence  it  is  accepted  or  rejected  according 
to  the  adopted  approval  criteria  [11], 

Other  convergence-enhancing  methods  can  also  be  used.  In  the  highly  non- 
linear cases,  the  initial  values  to  initiate  the  variable  vector  can  be  obtained  from 
a Poiscuillc  solution  which  can  be  easily  acquired  within  the  same  code.  The  con- 
vergence, as  indicated  already,  becomes  more  difficult  with  increasing  the  pressure 
gradient  across  the  flow  domain,  due  to  an  increase  in  the  nonlinearity.  An  effective 
approach  to  obtain  a solution  in  such  cases  is  to  step  up  through  a pressure  ladder 
by  gradual  increase  in  the  pressure  gradient  where  the  solution  obtained  from  one 
step  is  used  as  an  initial  value  for  the  next  step.  Although  this  usually  increases 
the  computational  cost,  the  increase  in  most  cases  is  not  substantial  because  the 
convergence  becomes  rapid  with  the  use  of  good  initial  values  that  are  close  to 
the  solution.  The  convergence  rate  and  speed  may  also  be  improved  by  adjusting 
the  flow  parameters.  Although  the  parameters  are  dependent  on  the  nature  of  the 
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physical  problem  and  hence  they  are  not  a matter  of  choice,  there  may  be  some 
freedom  in  tuning  some  non-critical  parameters.  In  particular,  adjusting  the  cor- 
rection factor  for  the  axial  momentum  flux,  a,  can  improve  the  convergence  and 
quality  of  solution.  The  rate  and  speed  of  convergence  may  also  depend  on  the 
employed  numerical  solver. 

Another  possible  convergence  trick  is  to  use  a large  error  margin  for  the  residual 
norm  to  obtain  an  approximate  solution  which  can  be  used  as  an  initial  guess  for 
a second  run  with  a smaller  error  margin.  On  repeating  this  process,  with  progres- 
sively reducing  the  error  margin,  a reasonably  accurate  solution  can  be  obtained 
eventually.  It  should  be  remarked  that  the  pore-scale  and  finite  element  models 
have  generally  different  convergence  behaviors  where  each  converges  better  than 
the  other  for  certain  flow  regimes  or  fluid-structure  physical  problems.  However, 
in  general  the  pore-scale  model  has  a better  convergence  behavior  with  a smaller 
error,  as  indicated  earlier.  These  issues,  however,  are  strongly  dependent  on  the 
implementation  and  practical  coding  aspects. 

6 Conclusions 

In  this  paper,  a pore-scale  modeling  method  has  been  proposed  and  used  to  obtain 
the  pressure  and  flow  fields  in  distensible  networks  and  porous  media.  This  method 
is  based  on  a residual  formulation  obtained  from  the  continuity  of  volumetric  flow 
rate  at  the  branching  junctions  with  a Newton- Raphson  iterative  numeric  technique 
for  solving  a system  of  simultaneous  non-linear  equations.  An  analytical  relation 
linking  the  flow  rate  in  distensible  tubes  to  the  boundary  pressures  is  exploited 
in  this  formulation.  This  flow  relation  is  based  on  a pressure-area  constitutive 
equation  derived  from  elastic  tube  deformability  characteristics. 

A comparison  between  the  proposed  pore-scale  network  modeling  approach  and 
the  traditional  one- dimensional  Navier-Stokes  finite  element  approach  has  also  been 
conducted  with  a main  conclusion  that  pore-scale  modeling  method  has  obvious 
practical  and  theoretical  advantages,  although  it  suffers  from  some  limitations  re- 
lated mainly  to  its  static  time  independent  nature.  We  therefore  believe  that 
although  the  pore-scale  modeling  approach  cannot  totally  replace  the  traditional 
methods  for  obtaining  the  flow  rate  and  pressure  fields  over  networks  of  intercon- 
nected deformable  ducts,  it  is  a valuable  addition  to  the  tools  used  in  such  flow 
simulation  studies. 
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Nomenclature 


a correction  factor  for  axial  momentum  flux 
f3  stiffness  coefficient  in  pressure-area  relation 
k viscosity  friction  coefficient 

H fluid  dynamic  viscosity 

v fluid  kinematic  viscosity 

p fluid  mass  density 

Poisson’s  ratio  of  tube  wall 

A tube  cross  sectional  area 

Ain  tube  cross  sectional  area  at  inlet 

Aa  tube  cross  sectional  area  at  reference  pressure 

Aou  tube  cross  sectional  area  at  outlet 

E Young’s  clastic  modulus  of  tube  wall 

/ flow  continuity  residual  function 

hQ  vessel  wall  thickness  at  reference  pressure 

J Jacobian  matrix 

L tube  length 

n number  of  network  nodes 

91  norm  of  residual  vector 

p pressure 

p pressure  vector 

Pi  inlet  pressure 

p0  outlet  pressure 

Ap  pressure  perturbation  vector 

Q volumetric  flow  rate 

r flow  continuity  residual 

r residual  vector 

R tube  radius 

t time 

x tube  axial  coordinate 
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